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Abstract 

The aim of our work is to study the effect of geometry variation on 
nucleation times and to address its role in the context of eukaryotic 
chemotaxis (i.e. the process which allows cells to identify and follow 
a gradient of chemical attractant). As a first step in this direction we 
study the nucleation dynamics of the 2d Ising model defined on a cylin- 
drical lattice whose radius changes as a function of time. Geometry 
variation is obtained by changing the relative value of the couplings 
between spins in the compactified (vertical) direction with respect to 
the horizontal one. This allows us to keep the lattice size unchanged 
and study in a single simulation the values of the compactification 
radius which change in time. We show, both with theoretical argu- 
ments and numerical simulations, that squeezing the geometry allows 
the system to speed up nucleation times even in presence of a very 
small energy gap between the stable and the metastable states. We 
then address the implications of our analysis for directional chemo- 
taxis. The initial steps of chemotaxis can be modelled as a nucleation 
process occurring on the cell membrane as a consequence of the ex- 
ternal chemical gradient (which plays the role of energy gap between 
the stable and metastable phases). In nature most of the cells modify 



their geometry by extending quasi-onedimensional protrusions (filopo- 
dia) so as to enhance their sensitivity to chemoattractant. Our results 
show that this geometry variation has indeed the effect of greatly de- 
creasing the timescale of the nucleation process even in presence of 
very small amounts of chemoattractants. 
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1 Introduction 



The decay of the metastable state at a first order phase transition is a rather 
well understood problem both in statistical mechanics and in quantum field 
theory. It attracted much interest in these last years due to its role in various 
different physical contexts ranging from the condensed matter applications 
to the study of the deconfinement transition in QCD or of the inflationary 
mechanism in cosmology The crucial problem in all these contexts is the 
evaluation of the lifetime of the metastable state from the knowledge of 
the microscopic interactions of the system and of the details of the non- 
equilibrium dynamics which drives it toward the true equilibrium state. If 
the Hamiltonian of the system only contains short range interactions then, 
as a general result, metastable states will always decay but their lifetime 
can vary by orders of magnitude. A simple but powerful effective model 
for the description of metastable decay is the so called droplet model which 
traces back to the seminal papers by Becker and Doring and Zeldovich [HE], 
and was later reinterpreted and extended using a field-theoretic language by 
Langer and Coleman and Callan [3H5] (for a review see for instance [6]). 
The main result of this approach is that the lifetime of the metastable state 
depends exponentially on the free energy F c of the so called "critical droplet" 
whose size R c can be evaluated by comparing the bulk and perimeter con- 
tributions to the free energy of the droplet. The weak metastability regime 
(in which we are interested in this paper) is reached when the F c is much 
larger than the temperature of the system (or than the Plank's constant in 
field theoretic applications). 

In this "classical" framework the system is assumed to be of infinite ex- 
tension and the background spacetime geometry is assumed to be flat and 
fixed. The picture becomes much more interesting if these two constraints 
are relaxed. When the system has a finite size L, comparable with the cri- 
tical radius R c , then the competition between these two scales may induce 
very different behaviours and creates a combination of different dynamical 
regimes [TJ. In particular for L << R c , the limiting scale of the nucleation 
process becomes the ratio L/£ where £ is the bulk correlation length of the 
system. When the background geometry is allowed to fluctuate then the 
exponential regime itself is abrogated and one finds instead a power law 
dependence of the decay rate [8] on the energy gap. 

The aim of our paper is to show, following the above observations, that 
by acting on the background geometry one can indeed modify the decay rate 
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of the system and reach, even for very small energy gaps between the stable 
and metastable phases, very short values for the lifetime of the metastable 
phase. 

Besides its intrinsic interest, this feature may be of some importance for 
biological applications. Nucleation theory, among several other biological 
applications, plays an important role in eukaryotic chemotaxis (the process 
which allows eukaryotic cells to identify and follow a gradient of chemical at- 
tract ant). Chemotaxis is mediated by two key enzymes: PI3K and PTEN. 
The standard way to model the process is to assume that the cell membrane, 
as a consequence of the external chemical gradient, undergoes a phase orde- 
ring process reaching a stable phase separation regime between a PI3K-rich 
phase (toward the chemoattractant) and a PTEN-iich phase (opposite to 
the chemoattractant) [9HTT]. In this picture the chemical gradient plays the 
role of energy gap between the stable and metastable phases. The limiting 
scale of this process is the cell size. Gradients smaller than a critical value 
would require critical droplets larger than the cell size to drive the process 
and this would lead, depending on the details of the model, to a random 
orientation of the two phases or to a dramatic increase in the timescale of the 
process. The way cells adopt to overcome this problem and enhance their sen- 
sitivity to chemoattractant is to modify their geometry by extending quasi- 
onedimensional protrusions (filopodia). This change of background geometry 
is exactly the type of process that we plan to study in this paper. 

All along the paper we shall use the two dimensional Ising model as 
prototypical example of a statistical system with short range interactions 
showing metastable phases when an external magnetic field is applied in the 
broken symmetry phase of the model. Besides this paradigmatic role there 
are two other important reasons which led us to concentrate in the Ising 
model. 

The first one is that it has been recently shown in [12] that the chemo- 
tacting model described above can be mapped into a suitable generalization 
of the 2d Ising model. The standard Ising model is obtained in the limit 
of infinite amount of cytosolic PTEN and PI3K enzymes (we shall discuss 
this mapping in some details below). 

The second reason is somehow more technical but of similar importance 
for our purposes. The 2d Ising model can be exactly solved (and the cor- 
relation length can be exactly evaluated) even for asymmetric values of the 
coupling constants. This fact will play a crucial role in the protocol that we 
shall use to modify the background geometry of the model. In fact it will 



2 



allow us to modify the geometry by simply acting on the coupling constants 
of the model without changing the structure of the 2d lattice (which would 
require a very sophisticated and inefficient software). 

This paper is organized as follows. In the next section we shall briefly 
discuss the chemotacting model mentioned above. In sect. 3 we shall intro- 
duce the Ising model with asymmetric couplings and briefly discuss its phase 
diagram. Sect. 4 will be devoted to a brief review of the classical nucleation 
theory (CNT) with a particular attention to the interplay between the two 
main scales of the problem: the system size L and the critical droplet radius 
R c . In sect. 5 we shall discuss the procedure we adopted to vary the geometry 
of the system and the results of our simulations while the last section will be 
devoted to a few concluding remarks. 

2 A simple chemotacting model and its Ising- 
like realization 

The link between nucleation dynamics in 2d Ising models and chemotaxis is 
based on the observation that chemotaxis requires as a preliminary step a 
phase separation process on the membrane of the cell [9HTT]. 

When a chemoattractant is switched-on in the external environment of 
the cell, owing to the interplay of the two enzymes PTEN (phosphatase and 
tensin-homolog) and PI3K (phosphatidylinositol 3- kinase), the phospholi- 
pids PIP2 (phosphatidylinositol bisphosphate) and PIP3 (phosphatidylinos- 
itol trisphosphate) are interconverted: PI3K catalyzes PIP2 phosphoryla- 
tion and PTEN catalyzes PIP3 dephosphorylation. While the phospholipids 
are bound to the inner surface of the membrane, the enzymes can freely dif- 
fuse in the cytoplasm, becoming active if absorbed on the membrane. When 
it happens, it is observed the formation of two patches in complementary 
domains, rich in PIP2 and PIP3 respectively. In particular, the PJP3-rich 
patch is localized on the side of the membrane exposed to the highest con- 
centration of chemoattractant (the leading edge) while the PIP2 rich-one on 
the other side (the rear edge). Similarly, a phase separation also involves 
cytoplasmatic localization of the enzymes PI3K and PTEN. 

Such spatial-organization phenomenon may be seen as a self-organized 
phase ordering process, where the cell state - driven by an external field - 
decays into two spatially localized chemical phases, thus defining a front and 
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a rear. This process lies at the heart of directional sensing, finally leading to 
cell movement toward the maximum of the attractant gradient. 

Following [12] the dynamics of the phase separation process occurring 
in the cell membrane can be modeled with a suitable generalization of the 
2d Ising model. The cell membrane is represented by a square lattice of 
side L with periodic boundary conditions and the spin variables Si — +1 
(— 1) are the enzymes PI3K (PTEN). The Hamiltonian of the model is 
composed by three terms: a short-range (i.e. involving nearest-neighborhood 
sites) ferromagnetic coupling describing the attractive interaction between 
enzymes, an external site-dependent magnetic field mapping the effect of the 
attractant and a long-range antiferromagnetic interaction term which keeps 
into account the finiteness of the cytosolic enzymatic reservoir. 

From the ratio of probabilities that PI3K and PTEN enzymes bind to 
the site z, it is possible to trace the energy difference between states S% — +1 
and Si — — 1, 



j€di 

(where j G di are the j nearest-neighborhood sites of the site i) which can be 
reinterpreted as the variation of the Hamiltonian 



where J is the ferromagnetic coupling constant, the first sum is over all the 
nearest-neighborhood sites (< i,j >), hi is the external site-dependent field 
representing the chemoattractant, N is the total number of enzymes PTEN 
and PI3K and A is a constant weighing the antiferromagnetic term with 
respect to the ferromagnetic one. The standard Ising model is recovered 
in the iV — > oo limit, i.e. for an infinite enzymatic reservoir. However it 
is important to stress that this is a non trivial limit. In fact due to the 
additional antiferromagnetic term the generalized 2d Ising model behaves 
in a quite different way with respect to the standard Ising model. In the 
broken symmetry phase it shows self-tuned phase separations also for very 
small amounts of gradient. However if the gradient is too small the two 
phases are randomly distributed on the lattice (i.e. the cell membrane) and 
chemotaxis cannot occur. In this generalized Ising model the absence of 
correlation between the spatial orientation of the two phases and the external 
field plays the role of the exponential increase of the mestable state lifetime 




(1) 



(2) 



<i,j> i i<j 
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in the standard Ising model. In this paper we decided to concentrate in the 
N — » oo Ising limit since it allows a rigorous theoretical description. We plan 
to devote a forthcoming paper to a detailed numerical analysis of its finite 
N generalization. 

In the above model the geometry plays no role: the lattice is a flat two 
dimensional torus which is assumed to be a reasonable approximation of a 
spherical cell membrane. However this is a very crude approximation from a 
biological point of view. Indeed in the chemotactic process, the geometry of 
the cell plays a pivotal role: signal sensing takes place as the leading edges 
of nearly one- dimensional finger-like projections (called filopodia) protrude 
from the cellular body. Such structures, as sophisticated antennas, allow the 
cell to sense even subtle concentrations of external signal, finally enabling 
motion (see [T3HT5] ). Since filopodia are dynamic structures appearing to 
sense the signal when the source of chemoattractant is far from the cell, the 
membrane is continually forced to change its geometry, squeezing and split- 
ting its leading edge. Moreover, since in a nearly one-dimensional geometry 
the interface tension between two different phases is lower than in a spherical 
geometry, even a very small gradient of external field can drive the coarsening 
toward its maximum. If the geometry fluctuates, the thin fingers developed 
from the leading edge can sense small amounts of gradient thanks to the 
rapid formation of clusters rich in PIP3 in the free extremities (and PIP2 in 
their complementary share) and then maintain the phase separation in the 
spherical geometry once they are reabsorbed by the membrane. 

As mentioned in the introduction the main goal of our work is to try 
to include these effects in the model by simulating the formation of a quasi 
onedimensional filopodium and then its reabsorption within the cell mem- 
brane. 

3 Ising model with asymmetric couplings 
3.1 General setting 

The 2d Ising model with asymmetric couplings is defined by the following 
Hamiltonian [IB] . 




(3) 



<i,j> 



<i,k> 



5 



where J x and J y denote the horizontal and vertical couplings respectively, 
the sums are over the horizontal (< i,j >) and vertical (< i,k >) nearest- 
neighborhood sites and the spins Si take as usual the values Si = ±1. We 
shall denote in the following with L x and 
horizontal and vertical directions respectively. 
The corresponding partition function is 



L y the sizes of the lattice in the 



^exp 

s 



<i,j> 



<i,k> 



(4) 



with K 



ks is the Boltzmann's constant and N 



L x L y 



k B T k B T ' 

is the number of sites of the lattice. 

The Kramer- Wannier duality holds also for this model and has the effect of 
relating the two couplings (for a review see for instance [TElfTT]). The duality 
relations are: 



sinh (2iT) sinh (2W) = 1 , sinh (2W*) sinh (2K) = 1 , (5) 

or equivalently: tanhK* = e~ 2W , tanh W* = e~ 2K . Denoting with ip (K, W) 
the free energy of the model we have: 

i) (K*, W*) = if; (K, W) + ~ In (sinh (2K) sinh (2W)) , (6) 

from eq.(E]) we can easily obtain the selfdual line which in this case is also 
the critical line of the model: 

sinh (2K) sinh (2W) = 1 . (7) 

For K = W this equation becomes the well known selfdual condition which 
allows to obtain the critical temperature of the standard 2d Ising model. 

All the points along this critical line correspond to the same universality 
class and are described in the continuum limit by the same 2d conformal field 
theory (CFT) with c = 1/2 (where c is conformal anomaly of the underlying 
CFT). 

Outside the critical line the equation sinh (2K) sinh (2W) = k, with k ^ 1 
constant, defines an infinite set of curves of constant temperature T in the 
K, W plane. In particular for k < 1 we have the T > T c (high-temperature) 
phase while for k > 1 we have T < T c (low-temperature) phase, where T c is 
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the critical temperature of the symmetric Ising model. In the following we 
shall be interested in the low T phase (i.e. k > 1) in which the Z 2 symmetry 
of the model is spontaneously broken. 

For reasons which will be clear below we choose different boundary con- 
ditions in the two directions: periodic b.c. in the vertical (y-axis) direction 
and free b.c. in the horizontal (x-axis) direction. The lattice thus acquires a 
cylindrical geometrjfj]. 

The most interesting feature of this model is that it can be solved exactly 
for any value of K and W in absence of an external magnetic field (see for 
instance [16]). Moreover an exact expression for the correlation length in 
the two spatial directions can be obtained in the high temperature phase, by 
exact resummation of the strong coupling expansion |17j : 

1 1 

6" = f T. . i „„u„.\i > iv = f TT^TTaT • ( 8 ) 
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tabnhK {l±taBhw 
i-tanhw 



In 



tanhW f i+tanhif 
\i-tanhii" 



From these expressions, using the duality transformation one can obtain 
the analogous values for the correlation length in the low T phase. This 
derivation however is not completely trivial. One should keep into account 
that in the vicinity of the critical point low and high T correlation lengths are 
related by an universal amplitude ratio = £+/£_ where £+ and £_ denote 
the correlation lengths above and below the critical point. In the case of the 
2d Ising universality class the value of is exactly known to be R^ = 2. This 
value keeps track of the fact that the spectrum of the underlying quantum 
field theory is very different in the two phases. 

These correlation lengths will play an important role in the following 
since they will allow us to evaluate the effective geometry of the model in 
the continuum limit. 



3.2 Asymmetric couplings and lattice deformations 

The most interesting feature of the Ising model with asymmetric couplings is 
that by changing the couplings along a line of constant absolute temperature, 

In the quantum field theory (QFT) interpretation of the model this choice of boundary 
conditions corresponds to the so called "finite temperature regularization" of the QFT. 
In this framework the model describes in the continuum limit a one dimensional QFT in 
contact with a heat bath whose temperature (which has nothing to do with the temperature 
which appears in the Hamiltonian of the model) is related to the inverse of the lattice size 
in the compactified direction (for a review see for instance |18|). 
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i.e. at fixed k, we can effectively modify the aspect ratio of the lattice. This 
is indeed a standard tool in finite temperature quantum field theory [T5]. In 
fact it is well known in this context that we may obtain two physically equi- 
valent regularizations of the model by squeezing the radius of the cylinder or, 
in a completely equivalent way, keeping the radius fixed and suitably increa- 
sing the coupling in the compactified (vertical) direction and simultaneously 
decreasing the coupling in the orthogonal (horizontal) directions. The exact 
amount of this change of the couplings can be determined using as physical 
scale the correlation length of the model. As a consequence of the changes in 
the couplings the vertical (i.e. along the compactified direction) correlation 
length t; y increases while the horizontal one £ x decreases. Since the only other 
scale of the problem is the circumference of the cyilinder L, models with the 
same ratio L/£ y are physically equivalent. We can thus trade a change in the 
lattice size in the compactified direction with a change in the value of the 
vertical coupling. This equivalence is routinely used in lattice gauge theories 
and has been recently confirmed also in the context of the 2d Ising model 
with a set of high precision simulations of the Binder cumulant [TH]. This 
equivalence will allow us to construct simulation protocols in which we can 
change in a continuous way the radius of the cylinder (keeping the lattice 
size fixed and changing the couplings) thus mimicking the elongation of eu- 
karyotic cells along the gradient of a chemoattracting potential and even the 
extreme situation of the formation of quasi one-dimensional filopodia. 



3.3 External magnetic field and metastability 



In order to study the metastable behaviour of the model we must add to 
the Hamiltonian an external magnetic field. In order to mimic an external 
chemotactic gradient, instead of choosing as usual an uniform external field, 
we choose to couple the system to a site dependent field whose shape is 
plotted in fig. flTJ. 

The Hamiltonian considered is therefore the following: 



H — J x ^ ] S{Sj J y *y ^ SiSk ^ ] hiSi 

<i,j> <i,k> i 



which leads to the partition function 



J^exp 

5 



K S l S J + W S t S k + pY,h l S l 

<i,j> <i,k> i 



(9) 



(10) 
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N, K and W as well as < i,j > and < i,k > are defined as above, (3 = -j^f 
and 



c is a constant, e defines the gradient intensity of the magnetic field and L x 
is the side of the lattice along the x-axis. With this definition the stable 
state is characterized by two phases (with minus sign for low values of x and 
plus sign for large values of x) separated by an interface while the metastable 
state, whose decay properties we shall study in the following, will be defined 
as the state in which all the spins are pointed in the —1 direction. 

4 Classical Nucleation Theory 

In this section we shall first briefly summarize the classical nucleation theory 
in the case of the Ising model defined on a square lattice of infinite size, 
then in the second part of the section we shall see the effect of introducing a 
finite size scale L comparable with the critical radius R c . Most of the results 
reported in this section can be found in standard textbooks. For a detailed 
discussion see [7]. 

4.1 Infinite systems 

In the Classical Nucleation Theory the decay of the metastable state is con- 
trolled by three main physical scales: the correlation length £ which sets the 
typical mean size of the droplet which are created in the system by ther- 
mal fluctuations, the critical droplet radius R c and the size R$ to which one 
droplet can grow before it is likely to meet another one. CNT usually as- 
sumes £ « R c « Rq and from simple thermodynamic arguments allows to 
obtain reliable estimates for the nucleation rate and for the mean lifetime of 
the metastable state. Let us briefly remind the main steps of this calculation. 

Let us first evaluate R c . To keep the analysis as simple as possible we 
shall assume the droplet to be of circular shape with area 7rR 2 and perimeter 




(11) 



is the external magnetic field, where 




(12) 
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2ixR. It is easy to show that the analysis holds for any generic shape. The 
free energy of the droplet is 

F(R) = 2nRa - nR 2 AE , (13) 

where AE is the difference in bulk free-energy density between the metastable 
and stable states and <7o is the surface tension. In the Ising case in which 
we are interested the energy gap is proportional to the absolute value of the 
magnetic field \H\. For not too large magnetic fields a good approxima- 
tion for AE turns out to be AE = 2MH/ (5 where M is the spontaneous 
magnetization. 

The probability of a fluctuation of this size is given by 

P rsj e . (14) 

Applying standard droplet theory arguments the critical radius is the 
value which minimizes (fT4l) 

Rc(T,H) = £L. (15) 
The nucleation rate (i.e. the probability of a critical fluctuation) is dominated 

2 

by the free energy of the critical droplet F c = ^| : 

I(T, H) ~ e~^r ~ e = e M (16) 

where we have introduced the constant T = -t^j 2 - which keeps track of 
all the microscopic details of the model. 

Let us discuss a few important features of this classical result. 

• Eq.( fl6|) is only the first order in the semiclassical approximation [3H5]. 
The next to leading order can be evaluated using quantum field theory 
methods and leads to a prefactor whose power can be predicted anali- 
tically and agrees very well with high precision Montecarlo simulations 
(see for instance [20]). We shall not further discuss this issue in the 
present paper since for our purposes the semiclassical result will be 
enough. 

• In eq.f ll6p the information about the first scale we have in the game, 
the correlation length of the model £, is hidden in the constant T. More 
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precisely in the scaling region the string tension is related to £ by the 
universal amplitude ratio R a = a £ [2T]. The requirement £ << R c is 
thus translated into the condition F » \H\ and allows us to have a 
precise estimate of the values of H for which the CNT can be used. At 
the same time the T » \H\ condition implies that we should expect 
exponentially long nucleation times. 

• In the above derivation we assumed a single droplet approximation. 
This is a reliable choice as far as critical droplets are far enough from 
each other. This is the physical meaning of the condition R c << R 
mentioned above. It can be shown (see for instance [7]) that Ro in- 
creases exponentially with l/|i?| thus, since R c ~ 1/\H\, there will 
always be a value of H below which this approximation is correct. 

4.2 Finite size systems 

If we consider a system of finite size L there will always be a value of \H\ 
small enough such that R c > L. In this limit the system behaves as along 
the coexistence line i.e. as if \H\ would be negligible with respect to the 
other scales. Thus we can neglect the area term in eq.f lT3|) and the droplet 
free energy becomes proportional to R. The saddle point analysis discussed 
above does not hold any more and one can show that instead the metastable 
lifetime increases exponentially with the system size L. To address this point 
let us define as the fraction of the system occupied by the droplet, i.e. 
(assuming again the simplified picture of a spherical droplet) <fi = ttR 2 /L 2 . 
Then we can rewrite the free energy as 



As above we can trade the line tension o"o for the correlation length using 
the universal amplitude ratio R a and rewrite F as 



with A = 2y/iuf)R a which, if we are interested in the free energy of a droplet 
which covers a finite fraction of the system, is a constant of order unity. Thus 
we see that the parameter which governs the droplet formation is actually the 
ratio L/£ and that the lifetime of the metastable state increases exponentially 
with 




(17) 



F(R) ~ 



(18) 
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In particular, if we are interested in a cylindrical geometry, the limiting 
scale will be the radius of the cilinder and, in case of asymmetric couplings, 
the correlation length in eq. ffTB"]) will be the correlation length along the 
compactified direction. 

5 Simulation setting and results 

Simulations were performed using a standard Glauber dynamics in order to 
mimic the local interactions of enzymes on the membrane surface |12j . 

We studied the model on a square lattice of sizes L x = L y = L = 100. 
The system was initialized setting all the spins to the S{ = —1 value. The 
boundary conditions were chosen to be periodic in the y-axis and free in the 
For defmiteness we set in eq. fllllll2p c = 1 and considered six values 
of gradient intensity, ranging from e = 5 to e = 0.0005 (see tabjl]). 

We simulated the system for three different values of k (see figj2]): k = 1.8 
(K = W = 0.55), k = 7 (K = W = 0.85) and k = 13.2 (K = W = 1). Recall 
that the relation between (K, W) and k is sinh (2K) sinh (2W) = k and k > 1 
corresponds to the broken symmetry phase. 

In order to study the effects of geometry variation on nucleation for each 
value of k we smoothly changed the values of K and W from K = W to 
the value of W and K corresponding to £ y — 50 = ^ (i.e. reaching a quasi 
one dimensional geometry), we let the system thermalize in this asymmetric 
point and then moved it back again to the symmetric K = W geometry. 

We checked that, for all the values of k that we studied, the results were 
essentially independent from the velocity of this transformation provided the 
thermalization time in the asymmetric points was large enough. Thus for 
simplicity we shall only report the results corresponding to two simulations 
protocols. 

The first one (protocol I) corresponds to a stepwise change in the value 
of the couplings. More precisely we started with 5 x 10 3 iterations with 
K = W in order to thermalize the system, we then changed in a stepwise 
way the couplings to the values corresponding to the maximal asymmetry 
(last column of tabfS]) and let the system thermalize for 5 x 10 4 iterations. 
Finally, again in a stepwise manner, we change the couplings back to the 
symmetric value and run the simulation for 10 4 iterations. 

The second one (protocol II) corresponds to a more gradual change of 
the couplings. In this case, after the first 5 x 10 3 iterations at K = W, the 
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system is gradually squeezed by simulating it for 10 4 iterations for each one 
of the three choices of asymmetric couplings reported in tabJTJ Then it is 
gradually moved back to the symmetric case following the same procedure 
and is finally allowed to thermalize again for 10 4 iterations at K — W . 

For all the values of e and k that we studied we then compared our 
results with those obtained by keeping the couplings (hence the geometry) 
unchanged (protocol III). We used this third set of simulations to identify 
the minimum value of e detectable by the system. Details on the parameters 
used in the simulations can be found in tabfSJ All the three protocols run for 
exactly 65 x 10 3 iterations and can thus be compared among them without 
bias. 

Due to the non trivial shape of the external magnetic field, besides the 
usual magnetization M — £V Si we used the following order parameter 
[12] to measure the polarization degree of the system: 

„= ^p^ (19) 
1 l2i \Ci-c\ 



that is 



a — — 



V ' (20) 



2 



ecos | ^ 



The results of our simulations are summarized in tabJT] and figs.( j3l4ll5"l) . 
Let us discuss our main findings. 

• In fig. ([3]) we report the behaviour of a for three values of k and e = 0.5 
as a function of the Montecarlo time, keeping K = W (protocol III). 
When K = W, for a fixed value of e the lifetime of the metastable 
state increases exponentially as a function of k. This is a well known 
result [6J and is in agreement with classical nucleation theory. In fact as 
A; increases the interfacial tension (and hence the free energy F c of the 
critical droplet) increases. For the smallest values of e and the largest 
values of k this time is much larger than our simulation time and as a 
consequence for these values of k and e we did not observe a decay of 
the metastable state (see tabJT]). In our biological interpretation these 
values of e are below the detectability threshold for the cell. 

• In figs.(jl]) and (J5J) we report the behaviour of a for k = 13.2 and 
several values of e as a function of the Montecarlo time for simulations 
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following protocols I and II respectively. By squeezing the lattice we see 
that the order parameter a for the three values of e reported gradually 
moves from ~ to ~ 0.5: i.e. at the end of the simulation the spin 
up cluster is exactly localized around the maximum of the magnetic 
field. We also report for comparison the result of the simulation with 
no change in the couplings. In this case due to the rather high value 
of k the metastable lifetime is exponentially long and in fact the order 
parameter keeps the value a = for the whole simulation time. By 
gradually decreasing e we could identify as e = 0.0005 the threshold 
at which even with the squeezing protocol the system could not detect 
the external magnetic gradient. This value should be compared with 
the analogous value e = 0.5 in absence of squeezing and shows that we 
have an enhancement of more than three orders of magnitude in the 
sensitivity of the system to external stimuli. 

The observed behaviour can be intuitively understood as follows. By 
changing the couplings in protocols I and II we change the value of the in- 
terfacial tension in the horizontal versus the vertical directions. This leads 
to the creation of asymmetric domains (vertical strips with our protocols) 
which play the role of nucleation drops. For large enough asymmetry these 
domains can wind around the vertical direction due to the periodic boundary 
conditions, thus leading to an effective dimensional reduction. In this regime 
the model behaves effectively as a one dimensional Ising model and the ver- 
tical strips start to aggregate as it happens for the spin kinks in the Id Ising 
model. This process is driven by the external field and can be very slow if the 
gradient e is small. However, even if it may not suffice to completely polarize 
the lattice (see the shape of the a parameter for the lowest e values in figs. 01]) 
and (jSJ)), it is enough to destabilize the original metastable configuration thus 
allowing the system, when it is moved back to the symmetric values of cou- 
plings, to reach the stable, polarized state. This is well described by fig.© 
in which snapshots of the spin configurations are plotted at different time 
during the simulation, both in the symmetric case and following protocol II. 

6 Concluding remarks 

Recent studies [9HTT] pointed out that chemotacting eukaryotic cells react to 
external chemical signals firstly polarizing at level of their inner membrane: 
the interplay between the enzymes PTEN and PI3K in the cytoplasm leads 
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to the formation of two complementary clusters of phospholipids (PIP2 and 
PIP3) on the membrane. If the concentration of chemoattractant is larger 
than a threshold value, then the P/P3-rich cluster is localized around the ma- 
ximum of the signal, determining internal cellular polarity. The above polari- 
zation is reached thanks to a phase separation process which can be cleverly 
mapped on a 2d-Ising-like model mimicking a spherical cell [T2]. However, 
dynamical geometry alterations play a pivotal role during the phase of signal 
caption. In particular it is well known that several types of eukaryotic cells 
use filopodia, which are quasi-one- dimensional finger-like actin protrusions 
exposed at the leading edge of migrating cells, as sensors of local external 
environment. As very subtle directional antennas, they are sofisticated tools 
which allow the cell to sense teeny quantity of signal otherwise negligible 
because below the threshold of detectability for a rigid spherical cell. 

In this paper, in order to describe the phase separation process occuring 
on the membrane of these filopodia we slightly modified the model proposed 
in [12] so as to be able to define it on a cylindrical lattice with a radius which 
may change as a function of time. In order to keep the exact solvability of the 
model we concentrated in the N — > 00 limit (see above) which corresponds 
to the standard anisotropic Ising model. 

The main assumption behind our approach is that the phase separation 
scenario proposed in [12] for the spherical cells holds also for filopodia, despite 
their extreme geometry. This assumption is supported by a few recent works 
on neuronal polarization showing that the P/3-kinase and its phospholipid 
product PIP3 are involved in the dendritic filopodial motility [22][23] . 

In particular clear evidences of a phase separation process were observed 
in [23], where the autors found that focal accumulation of PIP3 accompa- 
nies filopodial motility. In agreement with these findings we showed that 
squeezing the geometry of the cylinder (i.e. elongating the filopodia) allows 
to speed up the nucleation time even in presence of very small energy gaps 
between the stable and the metastable states. 

Altogether our analysis and these experimental results support the idea 
that eukaryotic cells use conformational changes of their membrane, and 
in particular the protrusion of filopodia, as a tool to optimize signal de- 
tection and speed up chemotaxis even in presence of very small amounts of 
chemoattractants by taking advantage of the enhanced efficiency of the phase 
separation process in a quasi-one-dimensional geometry. 
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Table 1: Summary of the system behaviour for each value of k and e. Legend: -k phase separation achieved 
with symmetrical couplings; • phase separation achieved with asymmetrical couplings; ® phase separation 
not achieved during the simulation time, either with asymmetrical couplings. 
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Table 2: Summary of the details on the parameters used in the simulations. For each value of k, the K and 
W values are reported together with the corresponding correlation lengths £ x and £ y . 




Figure 1: (Color online). Site-dependent external magnetic field for several values of e with c = 1. On the 
vertical axis is reported the site-dependent external magnetic field h(x) while on the the horizontal axis x 
runs on the L x lattice side. 
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Figure 2: (Color online). W = |asinh ^ sinh (2ic) y severa ^ vames °f On the V an d reported 

and K respectively. In the inset the points (K;W) corresponding to £ y = are plotted: the $ symbol 
refers to the point (0.0113, 2.53) for k — 1.8 , the # symbol refers to the point (0.0058, 3.55) for = 7 and 
the X symbol refers to the point (0.0054,3.9) for k = 13.2. 
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Figure 3: (Color online). Order parameters a and M as functions of Monte Carlo time for three values 
of k, e = 0.5 and symmetrical couplings K = W. From the top to the bottom of the figure a and M for 
k = 1.8 in black, for k = 7 in orange (grey) and for k = 13.2 in yellow (light grey). 
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Figure 4: (Color online). Order parameter a as function of Monte Carlo time for k = 13.2 and several values 
of e for a simulation following protocol I. From the top to the bottom of the figure, in orange (grey) data 
for e = 0.005, in black for e = 0.5, in red (dark grey) for e = 0.05 and in yellow (light grey) for e = 0.5 and 
symmetrical couplings. 
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Figure 5: (Color online). Order parameter a as function of Monte Carlo time for k = 13.2 and several 
values of e for a simulation following protocol II. From the top to the bottom of the figure, data for e = 0.5 
in black, for e = 0.05 in red (dark grey), for e = 0.005 in orange (grey) and for e = 0.5 and symmetrical 
couplings in yellow (light grey). 
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Figure 6: (Color online). The figure shows several snapshots of a square lattice with side L = 100 following 
protocol III (panel (a)) and protocol I (panel (b)). Each snapshot is taken at a particular Monte Carlo 
step, indicated by the number of the corresponding iteration (iter. no. 1, iter. no. 30000 and so on) below each 
square lattice. Blue (black) dots are up spins. In the panel (a), lattice with K — W for the whole simulation, 
k = 13.2 and e = 0.5. In the panel (b), coarsening dynamics in presence of asymmetrical couplings between 
5 x 10 3 and 55 x 10 3 iterations (from the second to the fifth lattice), k = 13.2 and e = 0.5. 



